#include "CHBDYGAREA3Calculator.h"
#include "Utility/Algebra/Tools.h"
#include "DataStructure/Input/Solver/ThermodynamicsConstant.h"

#include "DBManager/IService/Input/Node/INodeService.h"

#include "DBManager/IService/Input/Element/IElementService.h"
#include "DBManager/IService/Input/Subcase/ISubcaseService.h"
#include "DBManager/IService/Input/Material/IThermalMaterialService.h"
#include "DBManager/IService/Input/Property/IPCONVService.h"
#include "DBManager/IService/Input/Load/ICONVService.h"
#include "DBManager/IService/Input/Load/IQBDY3Service.h"
#include "DBManager/IService/Input/Constraint/ITEMPService.h"

#include "SGFem/Pre/Assembly/Assembly.h"
#include "SGFem/Pre/Assembly/ModelInfoQuery.h"

using SG::Algebra::Point;
using SG::Algebra::PointArray;
using SG::Algebra::Vector3D;
using SG::Algebra::Matrix;
using SG::Algebra::calcTriArea;
using SG::DataStructure::Common::Status;
using SG::DataStructure::Common::Real;
using SG::DataStructure::Common::Id;
using SG::DataStructure::FEM::ThermalMaterial;
using namespace SG::DataStructure::FEM;
using namespace SG::DBManager;
using SG::FEMSolver::AssemblyEigen;
using SG::FEMSolver::Element::CHBDYGAREA3Calculator;
std::size_t CHBDYGAREA3Calculator::s_nNode = 3;
std::size_t CHBDYGAREA3Calculator::s_nDim = 3;


void CHBDYGAREA3Calculator::Initialize (SG::DataStructure::FEM::ElementTypes type)
{
    m_type = type;
}

void CHBDYGAREA3Calculator::ComputeConvStiffness (const Matrix& coords, const ThermalMaterial& material,_OUT Matrix& ke2)
{
    Real h =  material.m_H;
    Real area;
    Vector3D normal;
    area = calcTriArea(coords, normal);

    ke2 = Matrix (3, 3, {2.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0});
    ke2 = ke2 * (h * area / 12.0);
}

void CHBDYGAREA3Calculator::ComputeConvStiffness (const SG::Algebra::Matrix& coords, const std::vector<Real>& h, _OUT SG::Algebra::Matrix& ke2)
{
    if (h.size() != 3) {
        // 丢出错误并终止程序
        throw std::runtime_error("Error: The Convetion boundary points of the area3 element should be 3.");
    }
    Real hsum =  (h[0] + h[1] + h[2]) / 3.0;
    Real area;
    Vector3D normal;
    area = calcTriArea(coords, normal);

    ke2 = Matrix (3, 3, {2.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0});
    ke2 = ke2 * (hsum * area / 12.0);
}

void CHBDYGAREA3Calculator::ComputeRadStiffness
    (
    const SG::Algebra::Matrix& coords, 
    const SG::DataStructure::Common::Real& Emissivity,
    const std::vector<SG::DataStructure::Common::Real>& T,
    const SG::DataStructure::Common::Real& Text,
    _OUT SG::Algebra::Matrix& ke4
    )
{
    if (T.size() != 3) {
        // 丢出错误并终止程序
        throw std::runtime_error("Error: The radiating boundary temperature points of the tetra element should be 3.");
    }
    
    std::vector<Real> hr(s_nNode);
    for (std::size_t i = 0; i < T.size(); ++i)
    {
        hr[i] = STEFAN_BOLTZMANN_NUMBER * Emissivity * (T[i] * T[i] + Text * Text) * (T[i] + Text);
    }

    ComputeConvStiffness(coords, hr, _OUT ke4);
}


void CHBDYGAREA3Calculator::ComputeConvLoad
    (
    const SG::Algebra::Matrix& coords, 
    const SG::DataStructure::FEM::ThermalMaterial& material,
    const SG::DataStructure::Common::Real& Text,
    _OUT SG::Algebra::Matrix& p3
    )
{
    Real h =  material.m_H;
    Real area;
    Vector3D normal;
    area = calcTriArea(coords, normal);

    p3 = Matrix (3, 1, {1.0, 1.0, 1.0});
    p3 = p3 * (h * Text * area / 3.0);
}

void CHBDYGAREA3Calculator::ComputeConvLoad
    (
    const SG::Algebra::Matrix& coords, 
    const std::vector<Real>& h,
    const SG::DataStructure::Common::Real& Text,
    _OUT SG::Algebra::Matrix& p3
    )
{
    if (h.size() != 3) {
        // 丢出错误并终止程序
        throw std::runtime_error("Error: The Convetion boundary points of the area3 element should be 3.");
    }
    Real hsum =  (h[0] + h[1] + h[2]) / 3.0;
    Real area;
    Vector3D normal;
    area = calcTriArea(coords, normal);

    p3 = Matrix (3, 1, {1.0, 1.0, 1.0});
    p3 = p3 * (hsum * Text * area / 3.0);

}

void CHBDYGAREA3Calculator::ComputeHeatFluxLoad (const Matrix& coords, const Real& q0, _OUT Matrix& p2)
{
    //Real qsum = (q0[0] + q0[1] +q0[2]) / 3.0;
    Real area;
    Vector3D normal;
    area = calcTriArea(coords, normal);

    p2 = Matrix (3, 1, {1.0, 1.0, 1.0});
    p2 = p2 * (q0 * area / 3.0);
}

void CHBDYGAREA3Calculator::ComputeElementMatries (const SG::DBManager::DBServiceFactorySP dbServiceFactory)
{
    auto& assembly = AssemblyEigen::GetInstance ();
    
    static const DofSet DOFs {DOFDirections::DOF_Ux};

    // SG::DataStructure::FEM::CTRIA3Data eleData;   //这个可能要用边界单元的导入
    PCONVData property; 
    ThermalMaterial material;
    CHBDYGData eleData;   
    // PCONVData property; 
    CONVData convData;
    TEMPData tempData;
    // ThermalMaterial material;
    std::vector<QBDY3Data> heatFluxDataArray;
    std::vector<TEMPData> tempDataArray;
    
    std::vector<std::size_t> nodeIndexes(s_nNode);
    std::vector<Id> dispCoordIds(s_nNode);
    std::vector<int> dofGlobalIndex(s_nNode);
    std::vector<Point> points(s_nNode);
    
    const auto& chbdygService = dbServiceFactory ->get<ICHBDYGService>();    //三角形边界单元数据格式
    const auto& nodeService = dbServiceFactory ->get<INodeService>();      
    const auto eleSize = chbdygService-> Size();
    for (std::size_t index{0}; index < eleSize; ++index)
    {
        /** 从数据库获取对应的单元数据 */
        chbdygService->FindByIndex (index, _OUT eleData);

        // Node coordinate
        ModelInfoQuery::QueryNode (nodeService, eleData.m_g, s_nNode, _OUT points, _OUT nodeIndexes, _OUT dispCoordIds);
        Matrix coords(s_nNode,s_nDim);
        for (std::size_t iNode = 0; iNode < s_nNode; ++iNode)
        {
            coords(iNode, 0) = points[iNode].m_x;
            coords(iNode, 1) = points[iNode].m_y;
            coords(iNode, 2) = points[iNode].m_z;
        }
        
        assembly.GetGlobalDOFId (nodeIndexes, DOFs, _OUT dofGlobalIndex);

        //获取材料参数   三角形边界单元实际上是四面体单元上的一个面，材料应该与四面体单元一致
        dbServiceFactory ->get<IPCONVService>()->FindById (eleData.m_id, _OUT property);
        dbServiceFactory ->get<IThermalMaterialService>()->FindById (property.m_MID, _OUT material);

        
        Matrix ke2(s_nNode, s_nNode);
        ComputeConvStiffness (coords, material, _OUT ke2);
        assembly.AccumulateStiffness (dofGlobalIndex, ke2);
        
        //dbServiceFactory ->get<IQBDY3Service>()   热流量q0
        Real q0 = 10.0;
        Matrix p2(s_nNode, 1);
        ComputeHeatFluxLoad (coords, q0, _OUT p2);
        //dbServiceFactory ->get<ICONVService>()    外部传热温度Text
        //dbServiceFactory ->get<ITEMPService>()
        // Real Text;
        Matrix p3(s_nNode, 1);
        // Id sid = 1.0;
        Id sid = 1;   //这个地方后期后改一改
        // ComputeConvLoad (coords, material, Text, _OUT p3);
        // assembly.AccumulateLoad (sid, dofGlobalIndex, p3 - p2);
   

        //获取对流换热系数和热交换温度
        if (Status::STATUS_Success == dbServiceFactory ->get<ICONVService>()->FindById (eleData.m_id, _OUT convData))
        {
            dbServiceFactory ->get<IPCONVService>()->FindById (convData.m_PCONID, _OUT property);
            dbServiceFactory ->get<IThermalMaterialService>()->FindById (property.m_MID, _OUT material);

            //对流刚度矩阵贡献
            Matrix ke2(s_nNode, s_nNode);
            ComputeConvStiffness (coords, material, _OUT ke2);
            assembly.AccumulateStiffness (dofGlobalIndex, ke2);
        
            //对流荷载贡献
            Matrix p3(s_nNode, 1);
            Real Text;
            dbServiceFactory ->get<ITEMPService>()->FindBySId(sid, _OUT tempDataArray);
            for (auto& tempData : tempDataArray)
            {
                if (tempData.m_G == convData.m_TA[0])
                {
                    Text = tempData.m_T;
                    break;
                }
            }
            ComputeConvLoad (coords, material, Text, _OUT p3);
            assembly.AccumulateLoad (sid, dofGlobalIndex, p3);
        }

        sid = 2;
        //获取热流参数
        if (Status::STATUS_Success == dbServiceFactory ->get<IQBDY3Service>()->FindBySId (sid, _OUT heatFluxDataArray))
        {   
            for (auto& heatFluxData : heatFluxDataArray)
            {
                if (eleData.m_id == heatFluxData.m_EID)
                {
                    //热流荷载贡献
                    Matrix p2(s_nNode, 1);
                    Real q0 = heatFluxData.m_Q0;
                    ComputeHeatFluxLoad(coords, q0, _OUT p2);
                    assembly.AccumulateLoad(sid, dofGlobalIndex, p2);
                }
            }
            continue;
        }

    }
}

// 后处理数据
void CHBDYGAREA3Calculator::DataRecovery (const SG::DBManager::DBServiceFactorySP dbServiceFactory, int subcaseId, SG::DataStructure::Common::Real time, const std::vector<SG::DataStructure::Common::Real>& globalDisplacement)
{
}

// 计算全局坐标系下单元应力刚度矩阵
void CHBDYGAREA3Calculator::ComputeStressStiffness (const SG::DBManager::DBServiceFactorySP dbServiceFactory, int staticSubId)
{
}
